#pragma once
#include <NR/nr3.h>
#include <NR/sort.h>
#include <NR/sparse.h>
#include <NR/linbcg.h>
#include <NR/asolve.h>

/* macros */
#define IX3(i,j,k) (((i)*(m_N+2)*(m_N+2))+((j)*(m_N+2))+k)
#define SWAP(x0,x) {float * tmp=x0;x0=x;x=tmp;}

class CFluidSolver
{
public:
	int m_N;
	float m_dt, m_diff, m_visc;
	int size;

	float *m_u, *m_v, *m_w, *m_u_prev, *m_v_prev, *m_w_prev;
	float *m_dens, *m_dens_prev;
	float *curl;

	//sparse matrix A
	NRsparseMat *matA;

public:
	CFluidSolver(int N, float dt, float diff, float visc);
	~CFluidSolver(void);

	void ClearData();
	void FreeMemory();

	//smoke simulation functions
	void LinearSolver(int b, float * x, float * x0, float a, float c);
	void PCBCGSolver(int b, float *x, float * x0, float a, float c);	//solving linear sys by preconditioned biconjugate gradient method

	void AddSource(float * x, float * x0);
	void Buoyancy(float * Fbuoy, float * temp);
	void VorticityConfinement(float * Fvc_x, float * Fvc_y, float * Fvc_z);
	void Diffuse(int b, float * c, float * c0, float diff);
	void Advect(int b, float * d, float * d0, float * du, float * dv, float * dw);
	
	//testing
	//void BFECC(int b, float * d, float * d0, float * du, float * dv);

	void Project(float * u, float * v, float * w, float * p, float * div);
	void SetBoundary(int b, float * x);

	void VelocitySolver();
	void DensitySolver();

	void SetSource(int i, int j, int k, float val);
};
